The imprints of primordial non-gaussianities on large-scale structure: scale dependent 

bias and abundance of virialized objects 
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We study the effect of primordial nongaussianity on large-scale structure, focusing upon the most 
massive virialized objects. Using analytic arguments and N-body simulations, we calculate the mass 
function and clustering of dark matter halos across a range of redshifts and levels of nongaussianity. 
We propose a simple fitting function for the mass function valid across the entire range of our 
simulations. We find pronounced effects of nongaussianity on the clustering of dark matter halos, 
leading to strongly scale-dependent bias. This suggests that the large-scale clustering of rare objects 
may provide a sensitive probe of primordial nongaussianity. We very roughly estimate that upcoming 
surveys can constrain nongaussianity at the level |/nl| 10, competitive with forecasted constraints 
from the microwave background. 



I. INTRODUCTION 



One of the fundamental predictions of standard (single- 
field, slow-roll) inflationary cosmology is that the density 
fluctuations in the early universe that seeded large-scale 
structure formation were nearly gaussian random (e.g. 
[l|, H, H, 0, Hi ) • Constraining or detecting non-gaussianity 
(NG) is therefore an important and basic test of the cos- 
mological model. To the extent that it can be measured, 
gaussianity has so far been confirmed; the tightest exist- 
ing constraints have been obtained from observations of 
the cosmic microwave background 0] . Recently, sev- 
eral inflationary models have been proposed which pre- 
dict a potentially observable level of nongaussianity, see 
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|22| . |23[ and [2J] for a review. Improved limits on NG 
would rule out some of these models; conversely, a ro- 
bust detection of primordial nongaussianity would dra- 
matically overturn standard inflationary cosmology and 
provide invaluable information about the nature of phys- 
ical processes in the early universe. In this regard, there 
has been a resurgence in studying increasingly more so- 
phisticated methods and algorithms to constrai n ( or, if 
we are lucky, detect) nongaussianity (25l. [26l. \2T\. \28L . 

Nongaussianity manifests itself not only in the cos- 
mic microwave background [s^, [HI, [13, [HI , but also in 
the late-time evolution of large-scale structure. For ex- 
ample, detailed measurements of higher order correla- 
tions like the bispectrum or trispectrum of galaxy clus- 
tering could provide a handle on primordial nongaus- 
sianity [H, [3^, [3^. The abundance of galaxy clus- 
ters, the largest virialized objects in the universe, has 
also long been recognized as a sensitive probe of pri- 
mordial NG [H, [13, M, M, \M,\M,\^- Because clus- 
ters are rare objects which form from the largest fluc- 
tuations on the tails of the density probability distribu- 
tion, their abundance is keenly sensitive to changes in the 



shape of the PDF such as those caused by nongaussian- 
ity. Large statistical samples of massive clusters have al- 
ready been compiled from wide-area optical imaging and 
spectroscopic surveys such as the Sloan Digital Sky Sur- 
vey [H, [4JI , the Two-Degree Survey [i^ , and from the 
Red Sequence Survey [11] and from X-ray surveys using 
the Chandra and XMM-Newton observatories [i^, [ll]. 
Future missions, such as the Dark Energy Survey, Su- 
pernova/Acceleration Probe and Large Synoptic Survey 
Telescope, will detect and study tens of thousands of clus- 
ters, revolutionizing our understanding of cluster physics 
as well as p roviding important constraints on cosmology 
[4i,[53,[5l|;[53,[5l,|53,[5l. 

To exploit the potential of these upcoming surveys as 
probes of primordial nongaussianity, it is important to 
calibrate the effects of NG on the abundance and cluster- 
ing of virialized objects. While no previous work has at- 
tempted to quantify the effects of NG on halo clustering, 
several groups over the past decade have constructed fit- 
ting formulae for the halo mass function [5^, [s^, [El] . All 
of this work, however, was a naly tic and relied on the va- 
lidity of the Press-Schechter [59j formalism, plus various 
further approximations. The resulting analytic estimates 
are, in general, rather cumbersome to compute and have 
questionable accuracy. As discussed below, the Press- 
Schechter model provides only a qualitative description 
of halo abundance, and fails to reproduce the halo mass 
function to within an order of magnitude over the mass 
and redshift range accessible to current and future clus- 
ter surveys. Therefore, analytic models for NG cluster 
abundance based on the Press-Schechter ansatz may not 
be sufficiently accurate. Given the high-quality data soon 
to be available, a much more precise calculation of clus- 
ter statistics will be required. Quite recently, two groups 
have attempted to quantify the mass function of clusters 
in NG models using N-body simulations [60,[6l|, reaching 
contradictory conclusions. 

In this paper, we use analytic arguments and numcri- 
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cal simulations to estimate the effect of NG on the abun- 
dance and clustering of virialized objects. Because N- 
body simulations can be expensive and there is a wide NG 
parameter space, we also strive to make our results use- 
ful to a cosmologist who is not necessarily equipped with 
the machinery or patience to run simulations or evaluate 
difficult analytic expressions. To this end. we provide a 
simple, physically motivated fitting formula for the halo 
mass function and halo bias, which we calibrate to our 
N-body simulations. 

Our main results are that the mass function and corre- 
lation function of massive halos can be significantly mod- 
ified by primordial nongaussianity. We find a somewhat 
weaker effect of NG on the mass function than previous 
analytic estimates. We also show analytically and numer- 
ically that NG strongly affects the clustering of rare ob- 
jects on large scales, implying that measurements of the 
large-scale power spectrum can place stringent bounds 
on NG. 

The plan of the paper is as follows. In Section [II] we 
derive analytic expressions for the abundance and clus- 
tering of rare peaks. In Section IIIII we describe our N- 
body simulations, followed in Section HVl bv a discussion 
of our measured halo mass function, and our fitting for- 
mula for the mass function. In Section |V] we present 
measurements of halo clustering within our simulations, 
and in Section IVll we discuss cosmological implications of 
our findings. 



II. ANALYTIC ESTIMATES 

In this section, we derive analytic expressions for the 
abundance and clustering of dark matter halos. As men- 
tioned above, such analytic approaches provide a use- 
ful qualitative framework for understanding gravitational 
collapse, however they cannot be used to describe quanti- 
tatively either the mass function or the clustering ampli- 
tude of collapsed objects. The expressions derived here 
are meant solely to motivate the more precise fitting for- 
mulae described in subsequent sections. 

We will focus on local NG of the form [J, [H, HI 



(1) 



In our notation, $ = — ^f, where is the usual New- 
tonian potential. On subhorizon scales, this choice of 
Newtonian gauge is valid, and the potentials $ and ^' 
satisfy the Poisson equation relating them to the over- 
density 5. On superhorizon scales, the Bardeen potential 
$ and ovcrdcnsity 5 are proportional, and not related by 
a Poisson equation, so our analysis will be valid only on 
subhorizon scales. With this choice of convention, posi- 
tive /nl corresponds to positive skewness of the density 
probability distribution, and hence an increased number 
of massive objects. 

For simplicity, we neglect the effect of the CDM trans- 
fer functions, which modify the shape of the $ power 
spectrum after nongaussianity is generated. Then the 



probability distribution for $ng is easy to write down, 
however the probability distribution for the density (Sng 
cannot be expressed analytically. Nevertheless, we can 
make progress by assuming that the NG correction is 
small, and by focusing only on high peaks of the density. 
The Laplacian of ^ng is 



V2$NG = V2(i 



2/nl[(/'V2 



(2) 



Because 0, V(/), and V^^ are all Gaussian fields whose 
statistics are fully specified by their power spectra, then 
Eqn. ^ above, relating (5ng = -(31^m/2ar|j)V^<I>NG 
to the Gaussian fields, allows us to determine fully the 
statistics of the nongaussian density (5ng- For example, 
the skewness of (5ng becomes, to lowest order in /nl, 
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On the average, the two terms 4'W'^<j> a-nd |V(/)p in 
Eqn. ([2]) are of the same order; the fact that they have 
equal but opposite expectation value is why (S-^q) = 
(6) = 0. However we are mainly interested in high peaks, 
where S oc — is large. Because |V0p is uncorre- 
lated with V^c/), and because at the peak of cj) its deriva- 
tive vanishes, we assume that |V^p may be neglected 
compared to ^V^^ in the vicinity of rare, high peaks. 
Then applying the Poisson equation near the peak gives 
(5ng ~ S[l + 2/nl0]- This expression appHes for the pri- 
mordial density and potential fields at early times. At 
late times, 5ng subsequently grows according to the lin- 
ear growth factor D{a), while the potential decays like 
g{a) oc D(a)/a. Therefore, rewriting this expression in 
terms of the late-time fields, we find 



(4) 



We see that the peak height is enhanced by a factor pro- 
portional to the primordial potential (f)p ~ <j>/g{a), rather 
than the evolved potential.^ 

Equation ([4]) will be the basis for the rest of our discus- 
sion. We emphasize that this is only valid in the vicinity 
of peaks, and so we focus on peaks for the remainder of 
this discussion. Because the fields S and (f> are Gaussian 
distributed, we can immediately derive properties of the 
distribution of Sj^q. For example, consider the mean shift 
in peak height for a peak of Gaussian density S: 



M 1 + 2/i 



NL 



(5) 



If the peak height S and background potential were 
uncorrelated, then there would be no systematic shift in 



^ An earlier version of this paper neglected to distinguish between 
the primordial and late-time potential, and hence omitted the 
g{a) factor. We are grateful to N. Afshordi for pointing this out 
to us. 
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peak height, and hence no change in the abundance of 
massive halos. However, S and 4> are correlated, imply- 
ing that rare peaks are systematically raised or lowered, 
depending upon the sign of /nl- Therefore, we expect 
changes in the mass function and the correlation func- 
tion. 

In the appendix, we derive expressions for the abun- 
dance and clustering of regions above a given threshold, 
which then give the clustering and mass function of halos 
in the Press-Schechter model. However, we can derive the 
form of the halo correlation function using a very simple 
argument. The halo correlation function is usually pa- 
rameterized in terms of the halo bias b, which is the rate 
of change of the halo abundance as the background den- 
sity is varied. Writing the matter overdensity as S and 
the halo overdensity as S^, wc can define the halo bias as 



bS. 



(6) 



It is normally assumed that b —^ const on large scales, 
but we will not make this assumption here. Consider a 
long- wavelength mode, providing a background density 
perturbation S and corresponding potential fluctuation 
(j). In the absence on nongaussianity, this perturbation 
raises subthreshold peaks above threshold, and thereby 
enhances the abundance of super-threshold peaks by b^S, 
where 6l is the usual (Gaussian) Lagrangian bias. For 
nonzero /nl, the long- wavelength mode also enhances the 
peak height by 2/NL'/'p<5pk, and we will focus on peaks 
near threshold, such that (5pk ~ dc- This provides an 
additional enhancement factor, giving a total 



(7) 



In Fourier space, the potential and density modes are 
related by = (30m/2ar|jfc^)(5, and so we see that the 
nongaussian bias acquires a correction 



Ab{k) = 25l/nl<5c 



(8) 



where again b^ refers to the usual Lagrangian bias for 
halos of this mass with Gaussian fluctuations. The total 
Lagrangian bias is then ^^(fc) = 6l -I- A6(fc). 

Since we have been working with the clustering of 
peaks in the initial density distribution, the above ex- 
pression for the bias applies only to the early-time, La- 
grangian bias. Translating these results to late-time, Eu- 
lerian bias is straightforward, however. The bias of Eule- 
rian halos is simply b = \+bL '■ the excess of halos in some 
Eulcrian volume with overdensity 5 is b5 = bL^ -\- 5. The 
first term corresponds to the excess of peaks in the initial 
Lagrangian volume, which are advected into the Eulcrian 
volume. The second term arises because an Eulcrian vol- 
ume with overdensity 5 has 5 times more mass than an 
average volume, and therefore 5 times more peaks. 

In summary, local NG generates a scale-dependent cor- 
rection to the bias of galaxies and halos, of the form 



A6(fc) =2(6-1)/nl<5c 
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FIG. 1: Slice through simulation outputs at 2: = gener- 
ated with the same Fourier phases but with /nl =—5000, 
—500, 0, -1-500, -1-5000 respectively from top to bottom. Each 
slice is 375 Mpc wide, and 80 Mpc high and deep. 
We can easily match by eye much of the large scale struc- 
ture; for example, an overdense region sits on the left, while 
an underdense region (void) falls on the right, in all panels. 
Note that for positive /nl, overdense regions are more evolved 
and produce more clusters than their Gaussian counterparts, 
while underdense regions are less evolved (e.(?. grid lines are 
still visible). For negative /nl, underdense regions are more 
evolved, producing deeper voids, while overdense regions are 
less evolved, as illustrated by the grid lines apparent in the 
left of the top panel. 



where 6 here now refers to the Eulcrian bias of the tracer 
population. In subsequent sections, we show that this 
simple expression, despite the underlying assumptions 
and approximations in its derivation, matches surpris- 
ingly well the halo clustering measured in our numerical 
simulations. 



III. NUMERICAL SIMULATIONS 

We numerically simulate the growth of structure in 
nongaussian cosmologies using the adaptive P'^M par- 
allel N-body code GRACOS^ [63, [63|. Non-gaussian ini- 
tial conditions were generated using the following pro- 
cedure. First, we generated a Gaussian random poten- 
tial field 0(x) using a power-law power spectrum with a 
scalar (density) index Ug — 0.96, and normalized so that 



2a g{a)r'jjk^ 



(9) 
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(78 = 0.76 [6| when multiplied by the matter transfer 
function. Following Refs. [3, [H, [631, then computed 
the nongaussian potential $ by adding a quadratic cor- 
rection in configuration space, 



$(a;) = 0(a;) + /NL(0'-(0')). 



(10) 



Wc then multiplied <f> by matter transfer functions in 
Fourier space for r2,„ = 0.24, SIa = 0.76, and computed 
particle displacements and velocities using the Zeldovich 
approximation [6^ . 

One immediate drawback to this approach is that, due 
to the strong Fourier mode coupling generated by the 
/nl term, our results may be affected by the absence 
of modes below the fundamental frequency or above the 
Nyquist frequency of our simulation volume. All N-body 
simulations can cover only a finite dynamic range, and 
therefore have zero power outside of their /c-space vol- 
umes. For Gaussian simulations, this is believed not to 
be a serious defect, because mode coupling is unimpor- 
tant on linear scales, and on nonlinear scales, the mode 
coupling generally transfers power to small scales. In our 
case, however, the /nl term couples all the modes sam- 
pled in our simulation to all the modes absent in our 
simulation. We have performed rudimentary estimates 
of the magnitude of this effect, by running simulations in 
which we high-pass or low-pass filter the /nl correction, 
and do not observe significant changes in the overall be- 
havior. Strictly speaking, however, it must be borne in 
mind that our results apply only for power spectra that 
are non-vanishing only over the finite range covered by 
our simulation volume. 

We have performed several simulations using both 
Gaussian and nongaussian initial conditions. For each 
Gaussian realization, we construct non-Gaussian realiza- 
tions using the same Fourier phases, with various /nl, 
e-g- /nl = ±500, ±50, and ±5. We ran simulations from 
a starting expansion factor a = 0.02 until the present 
time, a = 1, using 512'^ particles in a box of sidelength 
L = 8QQh^^ Mpc. For these parameters, each particle 
has a mass nip = 2.52 x 10^^ h~^MQ, so that clusters 
with masses exceeding M > 10^'^h~^ Mq are resolved 
with N > 400 particles. Since we are interested mainly in 
the masses and positions of cluster-sized halos, and not 
their internal structure, we have not used high force res- 
olution: we employ a Plummer softening length I of 0.2 
times the mean interparticle spacing. We have checked 
that using higher force resolution {I half as large) does 
not appreciably change the mass function. All simula- 
tions were performed at the Sunnyvale cluster at CITA; 
depending upon the value of /nl, the simulations com- 
pleted in 2-3 hours each on typically 8-10 nodes. As a 
consistency check, wc have also run a small number of 
1024^^ particle simulations with the same particle mass 
and force softening as above, but with twice the box size. 
These larger runs typically completed in 18-20 hours on 
64 nodes. In Figure [U we plot slices through our simula- 
tion volume at redshift z = 0, and the effects of varying 
/nl are readily apparent. Large positive /nl accelerates 
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FIG. 2: Mass functions measured from simulations with vari- 
ous /nl and identical phases (3 sets of initial conditions were 
used for each /nl). The top panel shows the mass function 
as well as the Gaussian fitting formula (dashed yellow line) 
from Warren et al. [g^. The bottom panel shows the ratio 
between the measured /nl = Gaussian mass functions and 
the respective non-Gaussian ones. 



the evolution of overdense regions and retards the evolu- 
tion of underdense regions, while large negative /nl has 
precisely the opposite effect. 



IV. THE HALO MASS FUNCTION 

We constructed late-time halo catalogues at redshifts 
z = 1, 0.5, and using the friends-of-friends group finder 
[131, with linking length b ~ 0.2. For Gaussian simu- 
lations, the halo mass function constructed this way has 
been extensively calibrated [1^, [6^ . Resulting mass func- 
tions are plotted in Figured! 



A. A new fitting formula 

Having measured the halo mass function, we next 
would like to construct a fitting function along the lines 
of those used for Gaussian simulations [i^, . As men- 
tioned above, previous techniques for estimating the non- 
gaussian mass function have been based upon the Press- 
Schechter [5§| ansatz. Given that the Press-Schechter 
mass function fails to match the halo mass function to 
within an order of magnitude over the mass and redshift 
ranges of interest to us [1^, and given the lack of any 
physical basis to the Press-Schechter ansatz [g^, [t^I , we 
have instead adopted an alternative approach which we 
describe next. 

We start by noting that the halo mass function dn/dM 
has been precisely calibrated for Gaussian cosmologies. 
Consider a Gaussian realization of the density field, 



which at late times evolves to produce halos with mass 
function dn/dMQ. As we slowly vary /nl away from zero, 
the structures forming at late times also slowly vary (c.f. 
Figure [1]), producing a different mass spectrum dn/dMf. 
If we vary /nl slowly enough, we can track the change 
in mass and position for individual halos: i.e., for each 
halo of mass Mq for /nl = 0, we can uniquely identify a 
corresponding halo of mass Mf for /nl 7^ 0, as long as 
|/nl| is sufficiently small. Since we know precisely the 
number of halos as a function of Mq, if we can determine 
the mapping Mq — > Mf, we will then have an estimate 
of the non-Gaussian mass function dn/dMf via 



dn 
dM} 



dMn 



dn dP 



dMn dM 



f 



{Mq), 



(11) 



where dP/dMf{MQ) is the probability distribution that a 
Gaussian halo of mass Mq maps to a non-Gaussian halo of 
mass Mf. Note that the probability distribution function 
dP/ dM f need not integrate to unity, / dAIf dP/ dM / 7^ 1 
in general, since the total number of halos is not con- 
served: halos can merge or split as /nl is varied. 

The next step is to determine the probability distribu- 
tion dP/dMf{Mo), by matching halos between Gaussian 
and non-Gaussian simulations. We match halos by re- 
quiring that matching pairs have significantly overlap- 
ping Lagrangian volumes; i.e. by requiring that halos 
have many particles in common, where particles are la- 
beled by their Lagrangian coordinates in the initial con- 
ditions. For each halo Mf in a non-Gaussian run, we 
loop over the halo's particles and identify which Gaus- 
sian halos own those particles in the run with /nl = 0. 
The Gaussian halo owning the largest fraction (exceed- 
ing 1/3) of the particles is then identified as the match 
for non-Gaussian halo Mf. Each Gaussian halo Mq can 
have one, several, or zero matching non-Gaussian halos, 
depending on /nl • By stacking Gaussian halos of similar 
mass Mq, we can determine dP/dMf{MQ). 

Examples of the probability distribution are shown in 
Figure SI and the mean and variance of the PDF are 
plotted in Figure [S] The behavior of the mean (Mf) and 
variance are quite regular, and appear consistent with 
simple power laws: 



Ml 

Mq 



1 



1.3 10-VNLfT8a(A/o,z)-2 (12) 
(S) " 1-4 10-'(/nl^8)"-'t(A/o,z)-\(13) 



where the rms ovcrdcnsity dispersion (j{M, z) is defined 
as usual by 



dh 



(14) 



where we use a top-hat window W{x) — 3ji{x)/x for 
R = (3M/47rp™)i/3, and P{k) and are the matter 
power spectrum and energy density respectively. 
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FIG. 3: Distribution of Mf as a function of Mq for one /nl = 
+500 simulation. The average shift towards higher masses is 
clearly visible. 



Because we desire a simple fitting formula, we assume 
that we can approximate the PDF as a normalized Gaus- 
sian whose mean and variance are given above, even 
though the PDF shape is quite clearly nongaussian (c.f. 
Figure |4]). As we show below, however, even this crude 
approximation is sufficient to achieve the ~ 10% precision 
in the halo mass function provided by standard fitting 
formulae for Gaussian simulations [7l|. Then Eq. (|lip . 
together with dP/dMf{MQ) which is assumed to be a 
Gaussian with the mean and variance given in Eqs. (jl2p 
and (Uni), fully specify our fitting function. Essentially, 
we have written the NG mass function as a convolution 
of the Gaussian mass function with a Gaussian kernel. 



B. Review of previous fitting formulae 

Before showing a comparison of the simulated mass 
function to our proposed fitting formula, we first describe 
alternative fitting formulae previously suggested in the 
literature: the Extended Press-Schechter (EPS) formal- 
ism , and the model of Matarrese et al. [5^, hereafter 
MVJ]. 



1, Extended Press-Schechter 

The EPS formalism generalizes the widely used Press- 
Schechter [5^ model, which posits that the fraction of 
mass in collapsed objects is equal to twice the fraction 
of the volume occupied by density peaks exceeding some 
critical ovcrdcnsity 5c. (The factor of two arises from the 
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FIG. 4: The probability distribution that a Gaussian halo of mass Mo maps to a non-Gaussian halo of mass Mf, i.e. 
dP/dMf{Mo). This plot can be understood as a (binned) slice through Fig. [3] Here we show the measured dP/dMf{Mo) 
(solid lines) and Gaussian fit (dashed line) for various /nl in the mass bin 1 < Mo / 10^* Mq < 3. The left panel corresponds 
to /nl > and the right panel corresponds to /nl < 0. Note that both the width and mean value of the PDF vary with /nl. 
The probability distribution is clearly poorly fit by a Gaussian, however as discussed in the text, it provides an adequate fit 
given the precision with which we can determine the halo mass function from N-body simulations. Whereas the high mass tail 
for /nl > (left panel) indicates that many /nl ~ halos will merge into more massive ones, the low mass tail for /nl < 
(right panel) accounts for the disruption of /nl ~ halos into lighter ones. 
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FIG. 5; Measured mean (left) and rms dispersion (right) of the mass shift, Mf/Mo — 1, as a function of a{Mo,z). Note that 
measurements at various redshift outputs {z — 0, 0.5, 1) have been combined in this plot, and that the /nl scaling has been 
divided out. The dashed lines shows our fits to these moments, c.f. Eqs. (|12|l and (|13|) . 



so-called 'cloud-in-cloud' problem [631 ■) Therefore the 
collapsed fraction becomes 



F{> M) 



So 



P{5; M)dS 



(15) 
(16) 



= 2 / PoHdiy 

where the probability distribution P{5] M) that the den- 
sity smoothed on mass scale M equals 5 is simply Pq{v), 
the Gaussian distribution with zero mean and unit vari- 
ance for V — 5c/cr, with a{M) given by Eq. (|14p . Note 



that i^(0) = 1 for hierarchical cosmologies where a{M) 
diverges as M — > 0; that is, all matter is assumed to be 
in viriaUzed objects of some mass. 

The differential mass function may then be readily 
computed: 



M- 



dn 
dM 



dF 



dM 



(17) 



For a Gaussian PDF, the mass M enters the right-hand 
side of this expression only via the lower bound of the 
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integral, Sc/(j{M), so we immediately obtain 



dn 



d\nM 



PS 



M a 



diner 



d\nM 



(18) 



This gives the well-known Press-Schcchter mass function. 

A class of fitting formulae based on this approach, and 
loosely called 'Extended Press-Schechter' (though appar- 
ently unrelated to the work of Refs. [zl, [ill) attempt to 
generalize this argument using nongaussian PDFs. The 
most obvious way of making this generalization, i.e. in- 
serting the nongaussian PDF P{6] M) into Eq. P?|) . faces 
several immediate difficulties, however. First, the Press- 
Schcchtcr factor of 2 is no longer valid; rather the cloud- 
in-cloud correction will depend upon the specific form of 
nongaussianity. Second, the shape of the PDF P{5\ M) 
now depends upon M and so we cannot simply replace 
the derivative of the integral in Eq. pT|) by the inte- 
grand. Lastly, and prosaically, the Press-Schechter mass 
function does not, in fact, fit the halo mass function in 
N-body simulations well, and so starting from PS is guar- 
anteed to fail in fitting the nongaussian mass function. 

The approach adopted by many previous workers (e.g. 
[t^ ) has instead been to assume that, although Press- 
Schechter cannot be used to derive the Gaussian mass 
function, it may be used to compute the departure of the 
mass function from its Gaussian value, i.e. 



nciM, z) 



J^F^g{>M) 
' Fg{>M) 



(19) 



dM 



where F is given by Eq. (|15p . In this approach, the non- 
gaussian mass function is computed by multiplying the 
Gaussian mass function (not Press-Schechter, but Jenk- 
ins et al. [ill or Warren et al. by the above ratio. 

The EPS prediction for the halo mass function is there- 
fore given by the derivatives of the PDF tails given in 
Eq. (|19p above. To implement this prescription, we com- 
pute the PDF tails directly from the initial conditions of 
our simulations: at each rcdshift wc arc interested in, we 
integrate the linearly evolved PDF to compute F^g{M) 
at masses ranging from about 10^^ to lO^^M©. Then we 
compute the mean value of F]^g{M) averaged over 10 
independent N-body simulations. Finally, we fit a cubic 
spline through the (computed mean of) Fj^Q{ln M) and 
differentiate with respect to mass. Evaluation of this for- 
mula becomes extremely difficult at high masses, simply 
because the statistics of peaks at these high masses be- 
comes too noisy. 



2. MVJ 

The MVJ [sst mass function is a further approximation 
to the EPS model described above. Instead of numeri- 
cally computing the PDF and its tails for each /nl, it is 
assumed that the ratio in Eq. (|19p may be determined 
from the skewness of the PDF. The expression for the 
mass function becomes 16111 
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where {6")m is nth moment of the density field evalu- 
ated on the characteristic mass scale M, and S^^m is the 
skewness on that mass scale. 

The advantage of the MVJ formula is that it does not 
require specification of the PDF of the density field — 
however, it does require knowledge of the skewness S^^m- 
In this work, we compute the moments of the density 
field directly from our simulations at the starting epoch 
a ~ 0.02, then scale them with the linear growth func- 
tion to the desired epoch. We can then evaluate the 
MVJ expression in Eq. (PT|) . Unlike the EPS approach 
described above, the formula does not become intractable 
at high masses. On the other hand, MVJ becomes oner- 
ously expensive to calculate at low levels of nongaus- 
sianity (e.g. |/nl| ^ 100), simply because the scatter 
in the measured skewness from run to run becomes com- 
parable to that generated by primordial nongaussianity. 
To see this, note that the scatter in {S^)m is roughly 
(73 ~ aljy^lB/N for N = Mho^/M samples. Approx- 
imating {5^)m ^ 6/nlO'mO'0, then for a 512^ grid and 
mass scale of 200 cells, and taking acf, = ix 10~^, requir- 
ing {S^)M/a3 > 5 translates into |/nl| ^ 100. 



C. Results and comparison to previous work 

Figure [S] shows the ratios 'ting/'t-G for /nl = 500 (left 
panel) and —500 (right panel). Simulation values are 
denoted with error bars, colored black (z = 0), blue {z ^ 
0.5) and red (z = 1). To compute the error bars in 
the ratios, and taking account of the fact that ting and 
TT-G measurements are correlated, we have adopted the 
larger error of the two alone (rather than adding them in 
quadrature), which is the error in riG (j^ng) for /nl > 
(/nl < 0). The solid lines denote our fits explained in 
Sec. mil Dashed lines refer to the EPS results, while 
the dotted lines represent the MVJ fitting function. The 
results clearly indicate that, while the EPS and the MVJ 
functions mutually agree"^, they both overestimate the 



^ The agreement between the EPS and MVJ is even better when 
an alternative expression is used in for /nl > 0, as pointed 
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FIG. 6: Ratios of the NG to Gaussian mass functions as a function of mass and at redshifts 2 = (black), z = 0.5 (blue) and 
z = \ (red). Points with error bars denote results from our simulations. Solid lines represent our fitting formula. Dashed and 
dotted lines denote the EPS and MVJ fitting functions respectively. Note that the EPS and MVJ agree mutually, but both 
significantly overestimate the effects of nongaussianity. (The discontinuity of EPS and MVJ fitting functions at M ~ 2-10^* Mq 
is due to transition from a smaller simulation box to the larger one.) 



effects of nongaussianity as found by our simulations, at 
a level typically ^ 100% although dependent upon mass 
and redshift. 

This result appears to disagree with the work of 
Kang et al. [g^I, who find a large discrepancy between 
EPS/MVJ and their simulations' mass function, in the 
sense that their simulations show a much larger effect of 
nongaussianity than predicted by the EPS type formal- 
ism. However, as noted by these authors, their simula- 
tions used a rather small number of particles 128'^) in 
a volume nearly 20 x smaller than ours, so it is unclear 
how well they probe the statistics of the rare objects of 
interest to us. In contrast. Grossi et al. [6l| have found 
very good agreement between the MVJ formula and their 
simulations' results. While our fitting function is in mild 
disagreement with the MVJ fitting formula, it is unclear 
whether our simulations are in disagreement with the 
simulations of Grossi et al. [6l1 |. Their simulations used a 
somewhat different cosmological model (higher erg) than 
ours, they have plotted cumulative rather than differen- 
tial mass functions, and of course the error bars in both 
their plots and ours are considerable. 

In summary, we conclude that our simple fitting func- 
tion appears consistent with the measured mass function 
from our simulations to within ^ 10% over the entire 
range of masses and redshifts that we consider. Since 
this is the level of precision that various N-body codes 
agree with each other in the mass function (Tlj , we have 



out Grossi et al. [6l| : see their Eq. (4). We have not used this 
correction in our Fig. |6] 



not attempted to achieve better agreement. EPS-like fit- 
ting formulae, such as the model of MVJ [1^, appear to 
overestimate the effects of nongaussianity. The level of 
discrepancy increases with increasing mass and redshift. 



V. HALO CLUSTERING 

Beyond one-point statistics like the halo mass function, 
N-body simulations also allow us to compute higher or- 
der statistics like the correlation function or its Fourier 
transform, the power spectrum. As shown in Sec. [H 
wc expect nongaussianity to produce pronounced effects 
on the halo power spectrum, specifically in the form of 
scale-dependent halo bias on large scales. This may seem 
somewhat surprising, due to very general arguments pre- 
viously given in the literature that galaxy bias is ex- 
pected to be independent of scale in the linear regime 
[TSI . [tgI . W\ . We can summarize the argument as follows. 
Suppose that the halo overdensity is some deterministic 
function of the local matter overdensity, 5h = F{5). On 
large scales, where \5\ <^ 1, we can Taylor expand this 
fimction, (5/i = a + 6(5 + . . .. Keeping only the lowest or- 
der terms and requiring that {5h) ~ then gives Sh = bd, 
which is linear deterministic bias. The key assumption in 
this argument was locality; i.e. that the halo abundance 
is determined entirely by the local matter density. N- 
body simulations with Gaussian initial conditions have 
confirmed that halo bias tends to a constant on large 
scales well in the linear regime. 

Once we allow for primordial nongaussianity. however, 
the above argument need not hold. For example, in this 
paper we have considered NG of the form f^h^^, and 
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FIG. 7: Cross-power spectra for various /nl. The upper panel 
displays Phs{k), measured in our simulations at 2 = 1 for ha- 
los of mass 1.6 x lO^M© < M < 3.2 x 10"Mq. The solid line 
corresponds to the theoretical prediction for Pss with a fitted 
bias feo=3.25. We see a strongly scale-dependant correction to 
the bias for /nl 7^ 0, increasing towards small k (large scales). 
The bottom panel displays the ratio /NL)/&(fc, /nl = 0). 
The errors are computed from the scatter amongst our simu- 
lations and within the bins. Triangles correspond to our large 
(1024^ particle) simulations whereas diamonds correspond to 
our smaller (512'' particle) simulations. The dotted lines cor- 
respond to our expression for the bias dependence on /nl 
defined in Eq. 



note that the gravitational potential is a nonlocal quan- 
tity. Hence the locality-based argument above does not 
apply for this form of nongaussianity, and our derived 
scale-dependence of the bias is not surprising. The spe- 
cific form we have derived is particular to the quadratic, 
local form of NG that we have assumed, however we ex- 
pect any NG that couples density modes with potential 
modes will in general lead to scale-dependent bias. On 
the other hand, nongaussianity of the form /nl<^^ does 
not lead to scale-dependent bias. 

In order to test our prediction for the scale dependence 
of bias, we have computed halo bias in our N-body simu- 
lations by taking the ratio of the matter power spectrum 
Pss and the halo-matter cross spectrum P^s = (<5^<5)- We 
have used the cross spectrum rather than the halo auto 
spectrum because the former should be less sensitive to 
shot noise from the small number of halos compared to 
DM particles. We have checked, however, that using the 
halo auto-spectra to compute bias gives consistent re- 
sults as the cross-spectra; i.e. we find no evidence for 
stochasticity. Examples of the various power spectra and 
resulting bias factors are plotted in figure Fig. [T] 

As can be seen, we numerically confirm the form of 
the predicted scale dependence. Because we focus on the 
statistics of rare objects, the errors on bias from individ- 
ual simulations plotted in Fig. [5] is large. We therefore 
attempt to improve the statistics on the comparison by 



FIG. 8: Ratio of the bias shift Ab measured from our simula- 
tions to that predicted by Eqn. Q, using 5c = 1.686. Biases 
were computed from cross-spectra measured on 28 simulations 
with 5 various /nl (-500, -100, 100, 500), 3 various redshifts 
(2 = 0, 0.5, 1) and 5 halo mass bins. Note that at higher 
k, nonlinear evolution also generates scale dependence in the 
bias [73. 



combining the bias measurements from multiple simula- 
tions. Figure [5] plots the average ratio between the bias 
measured in our simulations and our analytic prediction 
Eqn. ([9]), using Sr. ~ 1.686 as predicted from the spherical 
collapse model [79| . In computing the average plotted in 
this figure, we used a uniform weighting across the dif- 
ferent simulations, redshifts, and mass bins. Alternative 
weightings can shift the results by ~ 10%, so we conser- 
vatively estimate the systematic error in our comparison 
to be 20%. The agreement between our numerical sim- 
ulation results and our predicted bias scale-dependence, 
Eqn. ([S]), is excellent and perhaps surprising. Naively, 
we might expect a somewhat larger collapse threshold 6c 
to apply, considering the ellipsoidal rather than spherical 
nature of the collapse of halos in this mass range [tO] • 



VI. COSMOLOGICAL CONSEQUENCES 

Having derived fitting formulae for the abundance and 
clustering of halos in NG models, we now investigate how 
well upcoming surveys may constrain /nl, and whether 
NG could possibly affect the constraints derived on other 
cosmological parameters. We focus on galaxy cluster sur- 
veys and redshift surveys. Cluster surveys aim to con- 
strain cosmological parameters, in particular dark energy 
parameters, by exploiting the exponential sensitivity of 
the galaxy cluster abundance on cosmology. Similarly, a 
major goal for upcoming redshift surveys is to constrain 
dark energy by localizing baryonic acoustic oscillation 
(BAO) features in the galaxy power spectrum at mul- 
tiple redshifts. Examples of upcoming surveys include 
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the Atacama Cosmology Telescope^, South Pole Tele- 
scope^, Dark Energy Survey^, WiggleZ*", Planck^, Su- 
perNova/ Acceleration Probe^, and the Large Synoptic 
Survey Telescope^". 

Because primordial nongaussianity affects both the 
abundance and power spectra of massive halos, both of 
these types of surveys will be well-suited for constraining 
NG. On the other hand, potential NG could in principle 
degrade the expected constraints on dark energy param- 
eters, due to possible degeneracies. We use the Fisher 
matrix formalism to extract errors on seven cosmological 
parameters as well as /nl- Our estimates are only illus- 
trative; accurate forecasts for specific surveys will require 
a more sophisticated analysis. 



A. Constraints from P{k): 
Galaxy surveys, BAO and ISW 

We can (crudely) estimate constraints on parameters 
{pi\ derived from measurements of the power spectrum 
by assuming that bandpowers are measured with errors 
6P = (P + n^^)/-y/m, where P is the power in a band 
of width dk centered at wavenumber fc, n is the number 
density of galaxies, and m is the number of independent 
Fourier modes sampled by the survey, roughly given by 
TO = {2ti'^)~^V k'^dk [s^l. Then the Fisher matrix can be 
written as 

yfe„,.„ ^p^^pJ \ n) 27r2 

For simplicity, in Eq. (j23p we use the linear theory power 
spectrum and number density corresponding to z = 0.5, 
and assume an all-sky, volume-limited survey extending 
to z = 0.7. We integrate over wavenumbers between 
fcmin = 10^^/i/Mpc and fcmax = 0.1/i/Mpc. The results 
are insensitive to fc^ax but depend strongly on fcmin- We 
believe the fcmin used here is optimistic but reasonable. 
At high k (small scales), late-time nonlinear evolution 
can also generate scale-dependent bias [ill, however the 
redshift and scale dependence of this effect is quite dis- 
tinctive from NG and we ignore it here. 

We assume that the target galaxies have properties 
similar to luminous red galaxies (LRGs) [8l|, with co- 
moving number density n = 4 x 10~*(/i~^Mpc)~^ and 
bias 6o = 2. Equation then gives estimated errors 
on /nl of fT(/NL) ~ 7, which compares well with fore- 
casted constraints on nongaussianity for Planck. 



* http:/ /wwwphy.princcton.cdu/act/ 
^ http://spt.uchicago.edu 

'° http:/ /www.darkcncrgysurvcy.org 

^ http:/ /astronomy.swin.edu.au/wigglez/WigglcZ/Welcomc.html 

* http://www.rssd.esa.int/Planck 
^ http://snap.lbl.gov 

" http://www.lsst.org 



Unsurprisingly, we find little degeneracy between /nl 
and other cosmological parameters, given its distinctive 
effect on the shape of the power spectrum. Accord- 
ingly, there is little reason to believe that BAO deter- 
minations of dark energy parameters will be biased by 
nongaussianity, especially since the scale dependence is 
small over the wavenumbers of interest for the BAO wig- 
gles. To quantify this effect we determine the acoustic 
peak position by looking at extrema of the ratio of the 
power spectra with baryons and the power spectrum with 
zero baryons [s^ . When multiplying the matter power 
spectrum with baryons by our scale dependant bias, we 
find that /nl = 100 would shift the first BAO peak at 
fc ~ 0.07ft./Mpc by 0.4% at z = 1, and has a considerably 
smaller effect at the higher BAO peaks. The magnitude 
of this effect is comparable to the effect of non-linear 
corrections to the power spectrum [s^ [sj] , although the 
NG effect is primarily important on large scales while 
nonlinearities are most important on small scales. In 
principle, NG and nonlinearities could conspire to lead 
to a ~ 1 — 2% bias in the dark energy equation of state 
parameter w inferred from BAO observations [ssj . and 
so a careful joint analysis allowing both for NG and non- 
linear corrections will be required, which should not be 
difficult. 

Another probe of Ph& on large scales is the 
cross-correlation between cosmic microwave background 
(CMB) temperature anisotropics and large scale struc- 
ture, due to the integrated Sachs-Wolfe (ISW) effect 
[ssl . [sgI . [s?! . First detections of the ISW effect from cross- 
correlations of WMAP with various large scale surveys 
have been obtained with reported detections at the 2-4(t 
level [11, ii, Sill, [11, Sip, [11 . A combined analysis 
yields a ~ 5tT detection [96| . whereas a cosmic variance 
limited measurement would allow a ~ 7.5cr detection for 
the currently favored ACDM cosmology, and a somewhat 
more significant detection if the dark energy equation of 
state parameter is smaller [93, 13 • The cross-correlation 
between large-scale structure and CMB is directly pro- 
portional to a weighted projection of the scale-dependant 
bias. Since the z and fc dependence of our bias is very 
specific, we do not expect it to be severely degenerate 
with other parameters affecting the amplitude of the ISW 
effect (mostly w and fi™ for a flat universe). We can 
thus translate the ISW detection level into constraints 
on /nl- For the sake of simplicity, we assume that the 
ISW signal comes from z ~ 1 and is dominated by the 
angular multipole i ~ 20, corresponding to a wavenum- 
ber k ~ 6.66 X 10~'^/i/Mpc at z = 1. [oj. According to 
Eq. © the current 3(t [ha) detections of ISW translate 
into upper limits on |/nl| of 123 (61) (Icr) assuming a bias 
6o = 2, as appropriate for LRGs [93. A prospective 7.5 
a detection would translate into |/nl| % 38 (Ic). These 
estimates are clearly very crude, but are likely correct at 
the order of magnitude level. In comparison, the current 
limit from CMB bispcctrum measurements from WMAP 
give -54 < /nl < 114 (95% CL) @ whereas Planck is 
expected to constrain |/nl| < 10 (1 cr) [1^. 
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In summary, the large-scale galaxy power spectrum ap- 
pears capable of constraining local NG quite stringently 
for surveys reaching ~Gpc scales: |/nl| ^ 10. ISW 
or BAO observations could provide somewhat weaker 
bounds on /nl, though of course any constraints they 
can provide would be independent of the CMB bispec- 
trum and therefore worthwhile. Our estimates of fore- 
casted bounds on /nl were rather crude, but given the 
encouraging results, a more sophisticated treatment for 
specific survey parameters appears warranted. 



B. Constraints from cluster counts 

We next consider how well upcoming cluster surveys 
can constrain /nl by measurements of the cluster mass 
function dn/dM. Other forms of n onga ussianity may 
also be constrained by these surveys [lOOj, but we focus 
on the /nl form. For our fiducial survey parameters, we 
consider a fixed, redshift-independent lower mass limit of 
Miim = 2 X 1O"M0 and assume redshift bins of width 
Az = 0.1 uniformly distributed between z = 0.1 and 
z = Zmax < 2.0. We simultaneously vary seven cosmo- 
logical parameters besides /nl: A, the normalization of 
the primordial power spectrum at fcfid = 0.002ft, Mpc~^; 
physical matter and baryon densities flmh^ and flbh^, 
spectral index n^, the sum of the neutrino masses m^, 
the matter energy density today relative to critical ^Im, 
and the equation of state parameter of dark energy w. 
We assume no mass information (which would improve 
our parameter constraints) but also no systematic errors 
(which would degrade the constraints). We further as- 
sume 5000 square degrees on the sky, roughly consis- 
tent with expectations for the Dark Energy Survey or 
the South Pole Telescope. The fiducial survey has about 
7000 clusters (for erg = 0.76 cosmology) and about 23,000 
for as = 0.9). We use WMAP3 cosmological param- 
eters in determining error forecasts. The mass power 
spectrum A^(fc, a) = fc''P(fc, a)/(27r^) is written as 



A2(fc,a) 
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k 



k 



g\a)T\k) (24) 



where T{k) is the transfer function adopted from Eisen- 
stein and Hu [s^, and the growth function g{a) is com- 
puted exactly by integrating the well known seco nd o rder 
differential equation for growth (e.g. Eq. (1) in |l0l| ). 

Following the results of section IIV A| we assume that 
the mass function may be written as 
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where the nongaussian correction is computed using ei- 
ther our fitting formula, or EPS for comparison. For a 
given mass function, the total number of objects in a 
redshift interval of width Az and centered at z is 

f.z+Az/2 



iV(z, Az) = n., 
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FIG. 9: Forecasted errors on /nl from measurement of the 
cluster mass function, as a function of the maximum extent 
of the cluster survey, Zmax. Black solid and dashed line show 
the marginalized and unmarginalized error using our fitting 
formulae for riNG / , while the red lines show the errors us- 
ing the EPS formalism. Even though we use the average of 
the nongaussian PDF tail over 10 simulations, the EPS result 
becomes noisy at z > 1.5 due to poor sampling of the tails. 
Furthermore, it is clear that the EPS errors on /nl under- 
estimate those based on our simulations by up to a factor of 
three. 



where fi. 



is the total solid angle covered by the sur- 



vey, n(z, Mmin) is the comoving density of clusters more 
massive than M,nin, and dV/dfldz is the comoving vol- 
ume element. We assume flsurvoy = 5000 square degrees, 
roughly consistent with expectations for the Dark Energy 
Survey or the South Pole Telescope. The fiducial survey 
has about 7000 clusters (for erg = 0.76 cosmology) and 
about 23,000 (for erg = 0.9). 

Assuming Poisson statis tics, the Fisher information 
matrix reads [lo3, [lOl, [Toi 



^clus 
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dNk dN, 



Nk{zk,Az) dpt dpj 



(27) 



where pi are the 8 cosmological parameters including /nl, 
Nk is the number of clusters kth redshift bin, and the 
sum runs over the redshift bins extending to maximal 



redshift 



We assume no mass information (which 



Some attention needs to be paid when taking the derivative with 
respect to /nl i as it is especially with this parameter that the as- 
sumption that the likelihood function Gaussian may be violated, 
leading to results that are weaker or stronger than a full like- 
lihood calculation would reveal. We explore different values of 
d/NL and find convergence at IcJ/nlI $ 30. We also find, however, 
that the sensitivity to /nl (and also to the normalization A) is 
slightly higher when dfNL > {dhiA> 0) than when d/NL < 
(din A < 0); this is expected as excess of rare objects like galaxy 
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would improve our parameter constraints) but also no 
systematic errors (which would degrade the constraints). 

Lastly, we add a Planck prior on the parameter set, 
neglecting forecasted constraints on /ml expected from 
future CMB bispectrum measurements (since we arc in- 
terested in the sensitivity to NG of cluster counts alone) . 
The full Fisher matrix is given by 

^ ^ ^clus _^^CMB_ (28) 

Fig. [5] shows the result of our Fisher matrix estimate, 
i.e. the forecasted errors on /nl as a function of the max- 
imum extent of the cluster survey, z^ax- Black solid and 
dashed lines show the marginalized and unmarginalizcd 
error using our fitting formulae for nNc/'^G: while the 
red lines show the errors using the EPS formalism. The 
former errors are clearly well behaved, and asymptote at 
high Zmax as expected since cluster abundance rapidly 
vanishes. On the other hand EPS-produced /nl errors 
disagree with the simulations by up to a factor of three. 
Moreover, even though we use the average of the non- 
gaussian PDF tail over 10 simulations, the EPS result 
becomes noisy at z ^ 1.5 due to poor sampling of the 
tails. The magnitude of the discrepancy between EPS 
estimates and our simulations appears similar even for 
the higher erg = 0.9 model. As with our estimates from 
power spectrum constraints, we do not find significant 
degeneracies between /nl and other cosmological param- 
eters (correlation coefficient < 0.5). 

Since we find a weaker effect on cluster abundance than 
previous formulae like EPS or MVJ, this implies that con- 
straints found by Sefusatti et al. [72|, who performed a 
similar Fisher matrix estimate but with somewhat differ- 
ent assumptions, will be weaker once the NG sensitivity 
is calibrated off simulations. Direct quantitative compar- 
ison to two other relevant papers, Kang et al. (60j and 
Grossi et al. [6l|, is however difficult since these authors 
do not compute cosmological parameter error estimates. 



VII. DISCUSSION 

We have quantified the effects of primordial nongaus- 
sianity on the abundance and power spectra of massive 
halos. Our two principal results are as follows. 

First, we have provided a new fitting formula for the 
halo mass function. The formula is based on match- 
ing halos in Gaussian and non-Gaussian simulations: for 
/nl > the corresponding halos are more massive than 



clusters provides more cosmological leverage than their absence. 
This implies that true constraints on /nl will be slightly asym- 
metric around our fiducial value of zero; therefore, we make sure 
to take two-sided derivatives with cJ/nl = i30. For the EPS 
function, we take the two-sided derivative with cJ/nl = ±50 (re- 
call we ran sinulations with |/nl| = 5, 50 and 500), and check 
that the results are similar, if noisier, if d/NL = i5 is used. 



in the Gaussian case, and vice versa. The formula is con- 
sistent with the measured mass function from our simula- 
tions to within ^ 10% over the entire range of masses and 
redshifts that we consider. Being essentially a convolu- 
tion of the Gaussian mass function and a Gaussian kernel 
(Eqs. pip - (fT3|) ). the formula is also easy to use and does 
not require estimating the extreme tails of the nongaus- 
sian PDF of the density field. Our results also indicate 
that previous work based on Extended Press-Schechter 
type formulae overestimated the effects of nongaussian- 
ity on the abundance of halos by a factor of ~ 2 over the 
relevant mass scales. 

Secondly, we showed both analytically and numerically 
that nongaussianity (in the /nl model) leads to strong 
scale dependence of the bias of dark matter halos. We 
find remarkably good agreement between our analytic ex- 
pression and our numerical results. Measurement of the 
power spectrum of biased objects therefore provides a 
new avenue to detect and measure nongaussianity. While 
cluster counts can constrain NG at a level comparable to 
existing CMB constraints, |/nl| ^ 100, we found that fu- 
ture large-scale redshift surveys can potentially do much 
better, roughly |/nl| ^ 10. We do not find significant 
degeneracies between /nl and dark energy parameters in 
our Fisher matrix calculations, cither for mass function 
measurements or power spectrum measurements. More 
precise estimates will require considerably more sophis- 
ticated treatments than we have attempted in our illus- 
trative examples above. 

We close this paper by considering, in light of our 
findings, the optimal methods for constraining NG of 
the /nl form. Measurements of the power spectrum 
would appear the most promising; observations of high 
redshift, highly clustered objects on large scales would 
allow the strongest constraints on the scale-dependent 
bias signature of /nl- Fortunately, upcoming BAO sur- 
veys will likely provide the necessary observations of, 
e.g. luminous red galaxies (LRGs). Photometric sur- 
veys may also be useful in this regard. Since the ef- 
fects of NG are most pronounced on large scales, rather 
than small scales, precise spectroscopic redshifts may not 
be necessary. Photometric redshifts with errors of order 
Az « 0.03 have already been achieved for LRGs and 
for optically selec ted group s and clusters with prominent 
red sequences f46l . ll05l[l06l |. At z = 0.5, this corresponds 
to roughly 100 /I'^Mpc comoving, fairly small compared 
to the '^Gpc scales where NG becomes most important. 
Since photometric surveys can cover wider areas more 
deeply than spectroscopic surveys, they may turn out to 
provide tighter bounds. 

Besides their abundance and clustering, the internal 
properties of massive halos may also be sensitive to non- 
gaussianity. For instance, the concentrations and sub- 
structure content of massive halos have been found to 
depend upon primordial NG |l07( . Our simulations 
lacked sufficient force resolution to explore this in de- 
tail, but we note in passing that multiple groups find 
a tension between observations of massive lensing clus- 
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ters and theoretical pred ictions for Gaussian perturba- 
tions [loi, [Toi, [ii3, [ml . 

Another intriguing possibihty for probing primordial 
NG is to use statistics of the largest voids in the uni- 
verse. Just as the abundance and clustering of high den- 
sity peaks are affected by nongaussianity, so are the same 
properties for deep voids (albeit with an opposite sign, 
c.f. Fig. [1]). In a sense, because voids are not as nonlin- 
ear as overdense regions, their properties are more eas- 
ily related to the initial Lagrangian underdcnsitics whose 
statistics are straightforward to compute. Voids may be 
detected at high redshift as a deficit of Lyman-a forest 
absorption features in QSO spectra. The Sloan Digi- 
tal Sky Survey (SDSS) has already measured spectra for 
high redshift QSO's over a roughly ^8000 de g^ ar ea, cor- 
responding to a volume of ^ 30(Gpc//i)"^ |ll2l |. Each 
QSO spectrum typically probes ~ 400/i~^ Mpc, and the 
typical transverse separation between QSO sightlincs in 
SDSS is ~ 100ft." 1 Mpc, (P. McDonald, priv. comm.) so 
measurements of the clustering of ^ 10 Mpc-sized voids 
on ~ Gpc scales may already be feasible. 

Finally, we note that our conclusions are based on sim- 
ulations implementing a very specific type of local pri- 
mordial nongaussianity quantified by the /nl parameter. 
The validity of our conclusions in the context of other 
type of primordial nongaussianity is the subject of ongo- 
ing studies. 



derived for Gaussian statistics, and then show how they 
are modified by /nl nongaussianity. 



1. Review of Gaussian results 

We begin by identifying massive halos at late times 
with high peaks in the initial density distribution S{x) = 
6p{x)/p. Note that we work entirely in early-time, 
Lagrangian coordinates x in this section rather than 
late-time, Eulcrian coordinates. Earlier work (sol . [69l . 
IllSl . Ill4l | has shown that the abundance of peaks above 
threshold 6c ~ 1.686 reasonably describes (at the order- 
of-magnitude level) the statistics of halos forming at sub- 
sequent times. We briefly review some of these previous 
results, as the methods will be used in our analysis. 

Following the Press and Schechter [5^ ansatz, we 
smooth the density field and assume that density peaks 
with 6 > 5c produce halos. The density smoothed on 
scale R is given by 



Sifix) = (27r)-^ / d-'k SkW{kR)e 



(A.l) 



where W{x) is some smoothing window, e.g. top-hat or 
Gaussian. Assuming that the Fourier modes 6k arc Gaus- 
sian distributed with power spectrum P{k) = 
then 5fj{x) also has a Gaussian distribution, with vari- 
ance 
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dink 



k^Pjk) 
27r2 



W'^ikR). (A.2) 



The probability Pi for a given randomly selected region 
to exceed the threshold Sc is then simply the integral of 
the Gaussian probability distribution. 



Pi 



°° dP 



ierfc ( 
2 \^/2 



dv{2TT) 



(A.3) 



where v — S/ag, and similarly Vc — 5c/ us- 

The same power spectrum P{k) describing the density 
variance crgiR) also gives the matter correlation func- 
tion ^(ri2), and, following an elegant argument by Kaiser 
[ll3l | , can also be used to determine the correlation func- 
tion of rare peaks. Let us compute the probability P2 
that two randomly selected regions separated by distance 
ri2 ^ R are both above threshold. Again, this is simply 
an integral over the (joint) Gaussian distribution: 



APPENDIX: THE ABUNDANCE AND 
CLUSTERING OF HIGH PEAKS 

In this appendix we derive analytic expressions for the 
abundance and clustering of regions above the spherical 
collapse threshold 5c- We first review previous results 



P7 



d5i 



d5: 



exp I 



27r|S|i/2 



(A.4) 



where S = ((Si, (S2), the covariance matrix is given by 

(A.5) 
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and the matter correlation function ^(''12) is given by 
fc3p(fc), 



din fc- 



27r2 



-W'ikR)joikn2). (A.6) 



Rescaling the (5's by their variance, this becomes 



P2 
where 



dv; 



exp I 



27r|SP/2 



1 ?/' 

1 



(A.7) 



(A., 



and wc follow the notation of BBKS |114j in writing the 
normalized correlation function as '0(ri2) = £,{fi2)/c1; 
note that i/' < 1- 

To evaluate the integral in Eq. (|A.7[) . we change co- 
ordinates to variables that are uncorrelated. Using the 
Cholesky decomposition of the covariance matrix S ~ 
C^C, we write u = C^y for new variables y. Next, 
we rotate coordinates y = R • a; to bring the point 
[vi = Vc, vi = Vc) along the x\ axis. Then the integral 
becomes 



P2 



27r 



e[(l,0)-C'r-R-a;-i/c] 



e[(0,l)-C'r-R-a;-i/c] 



(A.9) 



where the Heavisidc function accounts for the two in- 
tegration bounds. Since ?/; < Ij we can write this as 



1 

2^ 
1 

2^ 



J —c{xi — Xc) 



(A.IO) 



where c = ^/{l + - V), and Xc = Vc^2l{\ + -0). 

We could easily evaluate the integral for /(xi) in terms 
of the error function, but the resulting integral over x\ 
would then not be analytic. However, we can derive an 
approximate solution in the limit I'c S> 1. For integrals 
of the form 



(A.ll) 



we can construct an asymptotic series by repeated partial 
integrations : 



1 - 



1 



Xq \ Xq J 

(A.12) 

In our case, /(xc) = and f'{xc) ~ 2c. Therefore, in the 
limit S> 1, we obtain 



P2 



1 

— f 

27r 

1 

— f 

2n 



2c 



J_e-''c/(i+'A) 



+ 2 
{1 - ^py/-^"" 



(A.13) 



Comparing this expression to the probability for a single 
peak to be above threshold then gives the peak-peak cor- 
relation function ^pk, which in the limit i^c 3> 1, ^ 1 
becomes 



1+Cpk - P2/P? 



,i.,^(l-l/(H-VO) 



1 + 1^2^ - 1 + 



(A.14) 



and therefore the (Lagrangian) bias 6| = ^pk/^ becomes 
6l « '^c/Sc- (A.15) 



2. Nongaussianity 

Our discussion so far has merely reviewed previous re- 
sults for Gaussian fluctuations; we now turn to nongaus- 
sian fluctuations. As noted above, we focus on NG of the 
form 



NG 



= + /] 



(A.16) 



We adopt the approximation of SjT] that the heights of 
rare peaks are modified by NG as (5ng ~ ^[1 + 2/nl0]i 
where in this appendix we adopt the notation that (j) 
refers to the primordial potential. At late times, (j) decays 
as the growth suppression factor g{a). 

Let us first consider the one-point distribution of peaks 
above threshold, (5ng > ^c- We express this as an integral 
over Gaussian variables and S; the integration bound 
then becomes (5ng = (5(1 + 2/nl0) > ^c, and using the 
fact that typically /nl|0| ^ 1, we have 6 > (5c(l— 2/nl'/'). 
The probability for (5ng to exceed threshold then is 



Pi 



dd 



<5c(1-2/nl0) 
00 



dfi 



dv 



exp 



27r|S|i/2 



2^|S|i/2 



(A.17) 



where pi = 4>/cf<p-, v — 5/as, Vc ~ 5c/(J5, and 77 = 
S/nlO'^i'c. We write the off-diagonal part of the normal- 
ized covariance matrix S as (/ii/) = r, where the cross- 
correlation coefficient r is not to be confused with the 
peak-peak separation ri2 appearing above and below. By 
a coordinate transformation, we can orient the integra- 
tion bound along a single axis. Changing coordinates 
from (/i, v) to {fi,v = 1/ + r]fJ-), and noting that the vari- 



ance of V is 
to obtain 



Pi 



= 1 + 2r]r + T] , we can integrate out fi 



/27r 



dx e 



-x^/2 



1 



Xr 



2°*lvl) 



(A.18) 



v^/^J\ + 2'qr + 'if. The effect of NG on the 



where x 

peak abundance is therefore simply to rcscalc the thresh- 
old density by a mass- and rcdshift-dcpendcnt factor. 
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Next, we turn to the peak-peak correlation function. 
As in the Gaussian case, we write the probabihty for two 
points both to be above threshold as 



P-2 ^ J ^ " 
x6(;^2 



exp I 



(27r)2|S|i/2 

W2 - l^c), 



Vc) 

(A.19) 



where u = ((Ui, /i2 , , ^^2 ) 7 and the notation is otherwise 
the same as above. We write the off-diagonal parts of the 
normalized covariance matrix as {V1V2) = ip, {^1^2) = 7i 
(i^iAii) = ^'^ and {^1^.2) = P- As above, we change vari- 
ables from j/j to Vi = Vi -\- rj^i to align the integration 
bounds along the density coordinate axes. This allows 
us to integrate out the two potential variables, leav- 
ing behind a 2-D integral. Rescaling the remaining two 
variables by their (identical) variance and again writing 
~ \/\ + '2,r\r + 77^, brings the integral to the form 



exp (— ice • S • cc) 
' 27r|S|i/2 



(A.20) 



where the off-diagonal component of S is {X1X2) = X 
given by 



X 



tp + 2r]l3 + rf-f 
1 + 2?77' + rp 



(A.21) 



The form of equation (|A.20p is identical to Eq. (|A.7p . 
with Vc ^ Xc and tjj x- So we can immediately write 
down the approximate solution. 



J_e-^c/(i+x)™-2(l 
2tt 



-x) 



3/2 



(i-x)i/^ 



(A.22) 



Comparing with the single-peak probability, we obtain 
the peak-peak correlation function, in the limit I'c 3> 1, 



X< 1 



1 + Cpk = P2/PI ~ 1 



(A.23) 



which to lowest order in ?/ becomes 



277(/?^ 
'277/3) 



(A.24) 



where 6l was the Lagrangian bias obtained for Gaus- 
sian peaks, c.f. Eq. (jA.lsp . Note that in going from the 
first line to the second, we neglect tt/j relative to (5 since 
Ci55('"i2)/o'55 is smaller than i^,^i{r\2)lo^^ by a factor scal- 
ing like {RjrYif' , where R is the smoothing scale of the 
peak and ryi is the peak-peak separation. 

The peak-peak correlation function is now no longer 
simply proportional to the matter correlation function, 
implying that the peak bias is not independent of scale. 
Fourier transforming this expression gives the peak power 
spectrum, 

fpk = h\{Pi,& H- 4/nl<5cP05) (A.25) 
which gives a scale-dependent change in the bias due to 
NG of 



A6(A:) = 26l/nl(5c 



26l/nl<5c 



£V5 



2ag rjjk'^ 



(A.26) 



where we have used the relation between the potential- 
density cross-spectrum and the matter power spectrum 
P<t,s = {'i^m/'20'gr^k-^)Pss, arising from the Poisson 
equation. The total Lagrangian bias is then 5l(/s) — 

bL + Ab{k). 
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